#include <stdlib.h>
#include <stdio.h>
#include "nbody.h"
#include "video.h"
#include <time.h>
#include <math.h>
#include <string.h>
#include <sys/time.h>
#include <mpe.h>

double *x, *y; 
double *fx, *fy, *vx, *vy, *m;
int    nbodies;
int    width;
int    height;
int    nsteps;
double g;
double interval;
double min_dist;

int    **grid;

/** Creo las variables globales para la asignacion de cuerpos*/
int cuerposPorProceso;
int inicioBucle;
int numeroProcesos;
int yo;
int resto;
int primero;


/** Creacion de las variables para el calculo del tiempo, no es necesario para el funcionamiento del programa solo para la toma de tiempos */
double tComT;
double tCalcT;
double tComSt;
double tComF;
double tCalcSt;
double tCalcF;
int start_comp, stop_comp;

int init_bodies(char *file) {
     FILE *f; 
     int i;
 
     f = fopen(file,"r");
	
     fscanf(f,"%d %lf %lf %d %lf %d %d\n",&nbodies,&g,&interval,&nsteps,&min_dist,&width,&height);
     x  = malloc(sizeof(double) * nbodies);
     y  = malloc(sizeof(double) * nbodies);
     fx = malloc(sizeof(double) * nbodies);
     fy = malloc(sizeof(double) * nbodies);
     vx = malloc(sizeof(double) * nbodies);
     vy = malloc(sizeof(double) * nbodies);
     m  = malloc(sizeof(double) * nbodies);

	/** Inicializacion de las variables que usamos para medir el tiempo*/
	tComT=0;
	tCalcT=0;
	tComSt=0;
	tComF=0;
	tCalcSt=0;
	tCalcF=0;
     for (i = 0; i < nbodies; i++) {
		fscanf(f,"%lf %lf %lf %lf %lf\n",&x[i],&y[i],&vx[i],&vy[i],&m[i]);
     }
     fclose(f);
     printf ("NBODY PROBLEM____________________________________________________________________\n");
     printf ("nbodies: %d\tg var: %0.2lf\tinterval: %0.2lfsec\tnsteps: %d\tmin_dist: %0.2lfp\twidth %d\theight: %d\n", nbodies,g,interval,nsteps,min_dist,width,height);
     return 1;
}

int finalize_nbodies(char *file) {
     FILE *f;
     int i;
     
     f = fopen(file,"w");
     fprintf(f,"%d %lf %lf %d %lf %d %d\n",nbodies,g,interval, nsteps,min_dist,width,height);
     for (i = 0; i < nbodies; i++) {
                fprintf(f,"%lf %lf %lf %lf %lf\n",x[i],y[i],vx[i],vy[i],m[i]);
     }
     fclose(f);
     free(x);
     free(y);
     free(fx);
     free(fy);
     free(vx);
     free(vy);
     free(m);
     return 1;
}

double factorial(double x) {
 double fac = 1;
 int i;
 for (i=2; i<=x; i++) fac *= i;
 return fac;
}

double mypow(double x, double exp) {
 double result = 1;
 int i;

 result = 1;
 for (i=0; i<exp; i++) result *= x;
 return result;
}

double mysin(double alpha) {
      double result = 0;
      int i;

      for (i=0; i < MAX_ACCURACY;i++) {
            result += (double)(mypow(-1,i) / factorial((double)(2*i+1))) * (double)mypow(alpha, 2*i+1);
      }
      return result;
}

double mycos(double alpha) {
      double result = 0;
      int i;

      for (i=0; i < MAX_ACCURACY;i++) {
            result += (double)(mypow(-1,i) / factorial((double)(2*i))) * (double)mypow(alpha, 2*i);
      }
      return result;
}


double myatan(double alpha) {
      double result = 0;
      int i;

      for (i=0; i < MAX_ACCURACY;i++) {
            result += (double)(mypow(-1,i) / (double)(2*i+1)) * (double)mypow(alpha, 2*i+1);
      }
      return result;
}
int calc_nbodies() {
     printf("lo sentimos pero este no es el modo secuencial \n");

     exit (1);

}

int calc_nbodiesP() {
     double *ax, *ay;
	/* Creacion de las variables locales auxiliares */	
	double *fauxx, *fauxy;
     double fax, fay;
     double dist, f, alpha, aux;
	/* Creacion de las variables locales auxiliares*/
	double *AUx, *AUy; 
	double *Avx, *Avy;
     int i,j;   
     	
	/* Toma la medida inicial de tiempo */
	MPE_Log_event(start_comp, 0, "start comp");
	tCalcSt=MPI_Wtime();

	/* Inicializacion de las variables locales en las que guardaremos las fuerzas calculadas por cada proceso*/
	fauxx = malloc(sizeof(double) * nbodies);
     	fauxy = malloc(sizeof(double) * nbodies);
	for (i = 0; i < nbodies; i++) {
                fauxx[i]=0;
		fauxy[i]=0;
	}

     ax = malloc(sizeof(double) * nbodies);
     ay = malloc(sizeof(double) * nbodies);

	/* Se calcula el numero de cuerpos que se le asignara a cada proceso */
	cuerposPorProceso = (int) nbodies / numeroProcesos;
	/* Se calcula los cuerpos que no pueden ser asignados a ningun proceso sin ocasionar desiquilibrio, resto
 		es el resultado del modulo de los cuerpos por los procesos que tenemos */
	resto = nbodies % numeroProcesos;

	/* CREO Q NO HACE FALTA  */
	if (resto != 0){
		cuerposPorProceso = (int) nbodies / numeroProcesos;
		resto = nbodies % numeroProcesos;
	}
   
	/* Si el resto es 0, el numero de cuerpos por proceso es el mismo en todos y no hay porblemas*/
	if(resto==0){
		/* El bucle comienza en el primer cuerpo que afecta a cada proceso q es igual al producto del id del proceso por el numero de cuerpos por proceso*/
		/* Asignado el inicio y fin de bucle el funcionamiento del mismo es similar al secuencial, variando que ahora las fuerzas se guardan en las variables auxiliares*/
     		inicioBucle = yo * cuerposPorProceso;
		for (i = inicioBucle; (i < inicioBucle + cuerposPorProceso) && (i< nbodies); i++) {
			fx[i] = fy[i] = 0;
     		}
		
     		for (i = inicioBucle; (i < inicioBucle + cuerposPorProceso) && (i<nbodies); i++) {
       			for (j = i + 1; j < nbodies; j++) {
          			dist = sqrt((x[j] - x[i]) * (x[j] - x[i]) + (y[j] - y[i]) * (y[j] - y[i]));
          			if (dist > min_dist ) {
               				f = g * m[i] * m[j] / (dist * dist);
               				if (f > MAX_F) f = MAX_F;
		
               				aux  = (y[j]-y[i]) / (x[j]-x[i]);

	       				alpha = atan(aux);
		
               				fax   = f * mycos(alpha);
	       				fay   = f * mysin(alpha);
	       				fauxx[i] += fax;
               				fauxy[i] += fay;
               				fauxx[j] -= fax; 
               				fauxy[j] -= fay;
          			}
       			}
     		}
		/* Toma del tiempo fin de calculo del primer bloque*/
		tCalcF=MPI_Wtime();
		MPE_Log_event(stop_comp, 0, "end comp");
		/* Suma al teiempo de calculo acumulado el nuevo tiempo gastado*/
		tCalcT+=(tCalcF-tCalcSt);

		/* toma del tiempo de inicio de comunicacion*/
		tComSt=MPI_Wtime();	
		/* Se envias las variables locales al proceso padre quien las suma y conserva*/
		MPI_Reduce(fauxx,fx,nbodies,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
		MPI_Reduce(fauxy,fy,nbodies,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);

		/* liberamos la memoria consumida por las variables locales de fuerza*/
		free(fauxx);
		free(fauxy);

		/* tengo todas las fuerzas, se las envio a los demas para q hagan sus calculos */
		MPI_Bcast(fx,nbodies, MPI_DOUBLE,0,MPI_COMM_WORLD);
		MPI_Bcast(fy,nbodies, MPI_DOUBLE,0,MPI_COMM_WORLD);
		/* Se toma el tiempo fin de comunicacion*/
		tComF=MPI_Wtime();
		/* Se acumula a la variable de tiempo total de comunicacion el gastado en las operaciones anteriores*/
		tComT+=(tComF-tComSt);

		/* Toma del tiepo inicial de calculo de este bloque */
		MPE_Log_event(start_comp, 0, "start comp");
		tCalcSt=MPI_Wtime();
		/* Reserva de memora e inicializacion de las variables locales auxiliares usadas para la posicion y la velocidad de los cuerpos*/
		AUx = malloc(sizeof(double) * nbodies);
		AUy = malloc(sizeof(double) * nbodies);
	     	Avx = malloc(sizeof(double) * nbodies);
     		Avy = malloc(sizeof(double) * nbodies);
		for (i = 0; i < nbodies; i++) {	
        	        AUx[i]=0;
			AUy[i]=0;
	                Avx[i]=0;
			Avy[i]=0;
		}
		/* Inicializacion del bucle dependiendelo del proceso en lel que estemo pues el inicio sera el id del proceso por el numero de ciuerpos asignado a cada proceso*/
		inicioBucle=yo * cuerposPorProceso;
		/* CAmbiado el inicio y el fin del bucle el funcionamiento es siilar al secuencial pero ahora lor resultdos se escriben en las
 			variables locales usando las globales para obtener los datos */
	     	for (i = inicioBucle; (i< inicioBucle + cuerposPorProceso) && (i<nbodies); i++) {
          		ax[i] = fx[i] / m[i];
          		ay[i] = fy[i] / m[i];
          		Avx[i] = vx[i] + ax[i] * interval;
          		Avy[i] = vy[i] + ay[i] * interval;
          		AUx[i] = x[i] + Avx[i] * interval;
          		AUy[i] = y[i] + Avy[i] * interval;
          
          		if (AUx[i] <= 0) {
             			AUx[i] = 2;
             			Avx[i] = - Avx[i];
          		}
          		if (AUx[i] >= width ) {
             			AUx[i] = width - 2;
             			Avx[i] = - Avx[i];
          		}
          		if (AUy[i] <= 0) {
             			AUy[i] = 2;
             			Avy[i] = - Avy[i];
          		}
          		if (AUy[i] >= height) {
             			AUy[i] = height - 2;
             			Avy[i] = - Avy[i];
          		}
 		}
		/* toma del tiempo fin de operacion y suma a la variabe acumulativa que lleva este recuento*/
		tCalcF=MPI_Wtime();
		MPE_Log_event(stop_comp, 0, "end comp");
		tCalcT+=(tCalcF-tCalcSt);
		/* toma del tiempo de inicio de comunicaciones*/
		tComSt=MPI_Wtime();
		/* Enviamos las variables locales al proceso raiz para que las sume y las caonserve*/
		MPI_Reduce(AUx,x,nbodies,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
		MPI_Reduce(AUy,y,nbodies,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);

		MPI_Reduce(Avx,vx,nbodies,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
		MPI_Reduce(Avy,vy,nbodies,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
		/* Se envian los datos actualizados a los demas procesos */
		MPI_Bcast(x,nbodies, MPI_DOUBLE,0,MPI_COMM_WORLD);
		MPI_Bcast(y,nbodies, MPI_DOUBLE,0,MPI_COMM_WORLD);

		MPI_Bcast(vx,nbodies, MPI_DOUBLE,0,MPI_COMM_WORLD);
		MPI_Bcast(vy,nbodies, MPI_DOUBLE,0,MPI_COMM_WORLD);
		
		/* Toma del tiempo fin de comunicacion y suma al tiempo de cpmunicacion acumulado */
		tComF=MPI_Wtime();
		tComT+=(tComF-tComSt);
	} /* fin del if (resto = 0) */


/* 
 * AHORA EL FUNCIONAMIENTO ES SIMILAR AL ANTERIOR POR LO QUE PARA NO REPETIR UNICAMENTE COMENTAREMOS LAS DIFERENCIAS
 *
 * En esta seccion nos encargaremos de lo que ocurre cuando la division del numero de cuerpos entre el numero de procesos no es entera
 * es decir que aparece un resto que sera el valor del resto de esta division.
 */
	MPE_Log_event(start_comp, 0, "start comp");
	tCalcSt=MPI_Wtime();
	if(resto!=0){
		/* El reparto de los cuerpos sobrantes se hace de uno en uno por proceso y por orden, de modo que los procesos con un id
 		 	menor al resto tendran un cuerpo mas que los que no cumplen esta condicion*/
		if(yo<resto){
			/* ahora se inicializa el bucle en el numero de proceso * (cuerpos pr proceso +1) puesto que todos los procesos anteriores
 				al actual tienen un cuerpo de mas*/
			inicioBucle = yo * (cuerposPorProceso+1);
			for (i = inicioBucle; (i < inicioBucle + (cuerposPorProceso+1)) && (i< nbodies); i++) {
				fx[i] = fy[i] = 0;
     			}
     			for (i = inicioBucle; (i < inicioBucle + (cuerposPorProceso+1)) && (i<nbodies); i++) {
       				for (j = i + 1; j < nbodies; j++) {
          				dist = sqrt((x[j] - x[i]) * (x[j] - x[i]) + (y[j] - y[i]) * (y[j] - y[i]));
          				if (dist > min_dist ) {
               					f = g * m[i] * m[j] / (dist * dist);
               					if (f > MAX_F) f = MAX_F;
		
               					aux  = (y[j]-y[i]) / (x[j]-x[i]);

	       					alpha = atan(aux);
		
               					fax   = f * mycos(alpha);
	       					fay   = f * mysin(alpha);
	       					fauxx[i] += fax;
               					fauxy[i] += fay;
               					fauxx[j] -= fax; 
               					fauxy[j] -= fay;
          				}
       				}
     			}
		}else{
			/* Si el id del proceso es mayor que el resto el numero de cuerpos sera el habitual pero se inicializara
 				en el valor habitual incrementado en el resto puesto que esos suerpos ya se habran calculado*/
			inicioBucle = yo * cuerposPorProceso + resto;
				for (i = inicioBucle; (i < inicioBucle + cuerposPorProceso) && (i< nbodies); i++) {
					fx[i] = fy[i] = 0;
     				}
     			for (i = inicioBucle; (i < inicioBucle + cuerposPorProceso) && (i<nbodies); i++) {
       				for (j = i + 1; j < nbodies; j++) {
          			dist = sqrt((x[j] - x[i]) * (x[j] - x[i]) + (y[j] - y[i]) * (y[j] - y[i]));
          				if (dist > min_dist ) {
               					f = g * m[i] * m[j] / (dist * dist);
               					if (f > MAX_F) f = MAX_F;
		
               					aux  = (y[j]-y[i]) / (x[j]-x[i]);

	       					alpha = atan(aux);
		
               					fax   = f * mycos(alpha);
	       					fay   = f * mysin(alpha);
	       					fauxx[i] += fax;
               					fauxy[i] += fay;
               					fauxx[j] -= fax; 
               					fauxy[j] -= fay;
					}
       				}
     			}
		}
		tCalcF=MPI_Wtime();
		MPE_Log_event(stop_comp, 0, "end comp");
		tCalcT+=(tCalcF-tCalcSt);
		tComSt=MPI_Wtime();
	
		MPI_Reduce(fauxx,fx,nbodies,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
		MPI_Reduce(fauxy,fy,nbodies,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);

		free(fauxx);
		free(fauxy);

		/* tengo todas las f, se las envio a los demas para q hagan sus calculos */
	
		MPI_Bcast(fx,nbodies, MPI_DOUBLE,0,MPI_COMM_WORLD);
		MPI_Bcast(fy,nbodies, MPI_DOUBLE,0,MPI_COMM_WORLD);

		tComF=MPI_Wtime();
		tComT+=(tComF-tComSt);

		MPE_Log_event(start_comp, 0, "start comp");
		tCalcSt=MPI_Wtime();
		AUx = malloc(sizeof(double) * nbodies);
		AUy = malloc(sizeof(double) * nbodies);
	     	Avx = malloc(sizeof(double) * nbodies);
     		Avy = malloc(sizeof(double) * nbodies);

		for (i = 0; i < nbodies; i++) {	
        	        AUx[i]=0;
			AUy[i]=0;
	                Avx[i]=0;
			Avy[i]=0;
		}
		/* Se vuelven a cumplir las condiciones para la inicializacion del bucle*/
		if(yo<resto){
			inicioBucle=yo * (cuerposPorProceso+1);

		     	for (i = inicioBucle; (i< inicioBucle + (cuerposPorProceso+1)) && (i<nbodies); i++) {
          			ax[i] = fx[i] / m[i];
          			ay[i] = fy[i] / m[i];
          			Avx[i] = vx[i] + ax[i] * interval;
          			Avy[i] = vy[i] + ay[i] * interval;
          			AUx[i] = x[i] + Avx[i] * interval;
          			AUy[i] = y[i] + Avy[i] * interval;
          
          			if (AUx[i] <= 0) {
             				AUx[i] = 2;
             				Avx[i] = - Avx[i];
          			}
          			if (AUx[i] >= width ) {
             				AUx[i] = width - 2;
             				Avx[i] = - Avx[i];
          			}
          			if (AUy[i] <= 0) {
             				AUy[i] = 2;
             				Avy[i] = - Avy[i];
          			}
          			if (AUy[i] >= height) {
             				AUy[i] = height - 2;
             				Avy[i] = - Avy[i];
          			}
     			}
		}else{
			inicioBucle=yo * cuerposPorProceso +resto;

			for (i = inicioBucle; (i< inicioBucle + cuerposPorProceso) && (i<nbodies); i++) {
          			ax[i] = fx[i] / m[i];
          			ay[i] = fy[i] / m[i];
          			Avx[i] = vx[i] + ax[i] * interval;
          			Avy[i] = vy[i] + ay[i] * interval;
          			AUx[i] = x[i] + Avx[i] * interval;
          			AUy[i] = y[i] + Avy[i] * interval;
          
          			if (AUx[i] <= 0) {
             				AUx[i] = 2;
             				Avx[i] = - Avx[i];
          			}
          			if (AUx[i] >= width ) {
             				AUx[i] = width - 2;
             				Avx[i] = - Avx[i];
          			}
          			if (AUy[i] <= 0) {
             				AUy[i] = 2;
             				Avy[i] = - Avy[i];
          			}
          			if (AUy[i] >= height) {
             				AUy[i] = height - 2;
             				Avy[i] = - Avy[i];
          			}
     			}
		}

		tCalcF=MPI_Wtime();
		MPE_Log_event(stop_comp, 0, "end comp");
		tCalcT+=(tCalcF-tCalcSt);
		tComSt=MPI_Wtime();

		MPI_Reduce(AUx,x,nbodies,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
		MPI_Reduce(AUy,y,nbodies,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);

		MPI_Reduce(Avx,vx,nbodies,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
		MPI_Reduce(Avy,vy,nbodies,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);

		MPI_Bcast(x,nbodies, MPI_DOUBLE,0,MPI_COMM_WORLD);
		MPI_Bcast(y,nbodies, MPI_DOUBLE,0,MPI_COMM_WORLD);

		MPI_Bcast(vx,nbodies, MPI_DOUBLE,0,MPI_COMM_WORLD);
		MPI_Bcast(vy,nbodies, MPI_DOUBLE,0,MPI_COMM_WORLD);

		tComF=MPI_Wtime();
		tComT+=(tComF-tComSt);
 
	}


	/* liberamos la memoria reservada para las variables locales*/

	free(ax);
	free(ay);

	free(AUx);
	free(AUy);
	free(Avx);
	free(Avy);

     	return 1;
}

/** metodo usado para el calculo de los cuerpos mediante el reparto por cuerpos */
int calc_nbodiesC() {
     	double *ax, *ay;
     	double fax, fay;
     	double dist, f, alpha, aux;
	/**  creacion de las variables locales que se necesitan para los almacenamientos en paralelo */
	double *fauxx, *fauxy;
	double *AUx, *AUy; 
	double *Avx, *Avy;
	int i,j;
	int numPares;
	int valor1;

	/** Inicializacion de las variables de tiempo y de MPE */
	MPE_Log_event(start_comp, 0, "start comp");
	tCalcSt=MPI_Wtime();
	/** Inicializacion de las variables auxiliares*/
	fauxx = malloc(sizeof(double) * nbodies);
     	fauxy = malloc(sizeof(double) * nbodies);
	for (i = 0; i < nbodies; i++) {
                fauxx[i]=0;
		fauxy[i]=0;
	}

     	ax = malloc(sizeof(double) * nbodies);
     	ay = malloc(sizeof(double) * nbodies);

	/** numero de cuerpos que no se pueden asignar emparejados a cada proceso */
	resto = nbodies % (numeroProcesos*2);
	
	/** numero de parejas de cuerpos que se le asignan a cada proceso */
	numPares=(int) nbodies/(numeroProcesos*2);
	primero=numPares*numeroProcesos;

	/** repartimos los cuerpos entre los procesos de la forma primero-ultimo */
	for (i = 0; i<numPares; i++) {
		/**  valor1 es el primer cuerpo que cogera mi proceso */
		valor1=numeroProcesos*i + yo;
		fx[valor1] = fy[valor1] = 0;
		fx[nbodies-valor1 -1] = fy[nbodies-valor1 -1] = 0;
     	}

	/** calcula las fuerzas desde el cuerpo actual hacia todos los demas, menos las anteriores que ya estan calculadas */
     	for (i = 0; i<numPares; i++){
		/**  valor1 es el primer cuerpo que cogera mi proceso */
		valor1=numeroProcesos*i + yo;
		for (j = valor1; j < nbodies; j++) {
		/** j es el primer cuerpo que coge mi proceso por el extremo izquierdo de la lista de cuerpos */
          		dist = sqrt((x[j] - x[valor1]) * (x[j] - x[valor1]) + (y[j] - y[valor1]) * (y[j] - y[valor1]));
          		if (dist > min_dist ) {
              			 f = g * m[valor1] * m[j] / (dist * dist);
               			if (f > MAX_F) f = MAX_F;
               			aux  = (y[j]-y[valor1]) / (x[j]-x[valor1]);
	       			alpha = atan(aux);
               			fax   = f * mycos(alpha);
	       			fay   = f * mysin(alpha);
				/* rellenamos los valores auxiliares de cada proceso con los resultados de los cuerpos*/
	       			fauxx[valor1] += fax;
               			fauxy[valor1] += fay;
               			fauxx[j] -= fax; 
               			fauxy[j] -= fay;
          		}
       		}

		for (j = (nbodies-valor1-1); j < nbodies; j++) {
		/** ahora cogeremos los cuerpos corresspondientes por el extremo de la izquierda */
			dist = sqrt((x[j] - x[(nbodies-valor1-1)]) * (x[j] - x[(nbodies-valor1-1)]) + (y[j] - y[(nbodies-valor1-1)]) * (y[j] - y[(nbodies-valor1-1)]));
          		if (dist > min_dist ) {
               			f = g * m[(nbodies-valor1-1)] * m[j] / (dist * dist);
               			if (f > MAX_F) f = MAX_F;
               			aux  = (y[j]-y[(nbodies-valor1-1)]) / (x[j]-x[(nbodies-valor1-1)]);
	       			alpha = atan(aux);
               			fax   = f * mycos(alpha);
	       			fay   = f * mysin(alpha);
				/* rellenamos los valores auxiliares de cada proceso con los resultados de los cuerpos*/
	       			fauxx[(nbodies-valor1-1)] += fax;
               			fauxy[(nbodies-valor1-1)] += fay;
               			fauxx[j] -= fax; 
               			fauxy[j] -= fay;
          		}
       		}
     	}/**  fin del calculo de las fuerzas para los cuerpos con parejas */


	/** Si hay cuerpos sin pareja debido al numero de procesos */
	if(resto!=0){
		/** El priemr cuerpo que le toca al proceso sera igual al priemr cuerpo sin pareja mas el numero de proceso */
		i=primero+yo;
		/** Si el cuerpo obtenido esta sin pareja tenog q cogerlo y calcularlo, esta sin pareja si esta dentro del rango de los cuerpos del resto */
		if(i<primero+resto){
			for (j = i; j < nbodies; j++) {
          			dist = sqrt((x[j] - x[i]) * (x[j] - x[i]) + (y[j] - y[i]) * (y[j] - y[i]));
	          		if (dist > min_dist ) {
        	       			f = g * m[i] * m[j] / (dist * dist);
               				if (f > MAX_F) f = MAX_F;
               				aux  = (y[j]-y[i]) / (x[j]-x[i]);
	       				alpha = atan(aux);
               				fax   = f * mycos(alpha);
	       				fay   = f * mysin(alpha);
					/* rellenamos los valores auxiliares de cada proceso con los resultados de los cuerpos*/
	       				fauxx[i] += fax;
               				fauxy[i] += fay;
               				fauxx[j] -= fax; 
               				fauxy[j] -= fay;
          			}
       			}
		}
		/** si el resto es mayor que el numero de procesos habra procesos con mas de 1 cuerpo en esta reasignacion*/
		if(resto>numeroProcesos){
			/** cogo el cuerpo que este en la posicion del primer cuerpo sin pareja mas el numero de procesos mas el numero de mi proceso, solo si esta en el rango de cuerpos sin pareja */
			if((primero+yo+numeroProcesos)<(resto+primero)){
				i=primero+yo+numeroProcesos;
				for (j = i; j < nbodies; j++) {
          				dist = sqrt((x[j] - x[i]) * (x[j] - x[i]) + (y[j] - y[i]) * (y[j] - y[i]));
          				if (dist > min_dist ) {
               					f = g * m[i] * m[j] / (dist * dist);
	               				if (f > MAX_F) f = MAX_F;
        		       			aux  = (y[j]-y[i]) / (x[j]-x[i]);
	       					alpha = atan(aux);
	               				fax   = f * mycos(alpha);
		       				fay   = f * mysin(alpha);
						/* rellenamos los valores auxiliares de cada proceso con los resultados de los cuerpos*/
	       					fauxx[i] += fax;
               					fauxy[i] += fay;
               					fauxx[j] -= fax; 
	               				fauxy[j] -= fay;
        	  			}
       				}	
			}
		}
	}


	/** toma del tiempo fin de calculo, su acumulacion en la variable total y marca del evento para el jumpshot*/
	tCalcF=MPI_Wtime();
	MPE_Log_event(stop_comp, 0, "end comp");
	tCalcT+=(tCalcF-tCalcSt);

	/** Marca del comienzo del tiempo de comunicaciones*/
	tComSt=MPI_Wtime();
	/** envio de las variables locales al proceso raiz para el calculo global */
	MPI_Reduce(fauxx,fx,nbodies,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
	MPI_Reduce(fauxy,fy,nbodies,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);

	/** Se libera la memoria de las variables locales*/
	free(fauxx);
	free(fauxy);

	/* tengo todas las f, se las envio a los demas para q hagan sus calculos en un envio de Bcast*/
	MPI_Bcast(fx,nbodies, MPI_DOUBLE,0,MPI_COMM_WORLD);
	MPI_Bcast(fy,nbodies, MPI_DOUBLE,0,MPI_COMM_WORLD);
	
	/** Marca del fin del tiempo de las comunicaciones y acumulacion de este en la variable total*/
	tComF=MPI_Wtime();
	tComT+=(tComF-tComSt);

	/** Marca de inicio del tiempo de calculo*/
	MPE_Log_event(start_comp, 0, "start comp");
	tCalcSt=MPI_Wtime();

	/** reserva de la memoria para las variables locales auxiliares e inicializacion */
	AUx = malloc(sizeof(double) * nbodies);
	AUy = malloc(sizeof(double) * nbodies);
     	Avx = malloc(sizeof(double) * nbodies);
     	Avy = malloc(sizeof(double) * nbodies);
	for (i = 0; i < nbodies; i++) {	
                AUx[i]=0;
		AUy[i]=0;
                Avx[i]=0;
		Avy[i]=0;
	}

	/** recorro los cuerpos por el numero de parejas asignadas a cada proceso*/
	for (i = 0; i<numPares; i++){
		/** en cada iteracion cogo un elemento del inicio*/
		valor1=numeroProcesos*i + yo;

          	ax[valor1] = fx[valor1] / m[valor1];
          	ay[valor1] = fy[valor1] / m[valor1];
	        Avx[valor1] = vx[valor1] + ax[valor1] * interval;
	        Avy[valor1] = vy[valor1] + ay[valor1] * interval;
	        AUx[valor1] = x[valor1] + Avx[valor1] * interval;
        	AUy[valor1] = y[valor1] + Avy[valor1] * interval;
      
	        if (AUx[valor1] <= 0) {
        		AUx[valor1] = 2;
	             	Avx[valor1] = - Avx[valor1];
        	}
	        if (AUx[valor1] >= width ) {
        		AUx[valor1] = width - 2;
	             	Avx[valor1] = - Avx[valor1];
        	}
	        if (AUy[valor1] <= 0) {
        	     	AUy[valor1] = 2;
	             	Avy[valor1] = - Avy[valor1];
        	}
	        if (AUy[valor1] >= height) {
        	     	AUy[valor1] = height - 2;
	        	Avy[valor1] = - Avy[valor1];
        	}
		/** ahora cogo el elemento opuesto en la cola de cuerpos, es decir el q ocupa el mismo puesto en sentido inverso, siempre q sean cuerpos emparejados*/
		valor1=(nbodies - (numeroProcesos*i + yo))-1;

          	ax[valor1] = fx[valor1] / m[valor1];
          	ay[valor1] = fy[valor1] / m[valor1];
	        Avx[valor1] = vx[valor1] + ax[valor1] * interval;
	        Avy[valor1] = vy[valor1] + ay[valor1] * interval;
	        AUx[valor1] = x[valor1] + Avx[valor1] * interval;
        	AUy[valor1] = y[valor1] + Avy[valor1] * interval;
          
	        if (AUx[valor1] <= 0) {
        	     	AUx[valor1] = 2;
	             	Avx[valor1] = - Avx[valor1];
        	}
	        if (AUx[valor1] >= width ) {
        	     	AUx[valor1] = width - 2;
	             	Avx[valor1] = - Avx[valor1];
        	}
	        if (AUy[valor1] <= 0) {
        	     	AUy[valor1] = 2;
	             	Avy[valor1] = - Avy[valor1];
        	}
	        if (AUy[valor1] >= height) {
        	     	AUy[valor1] = height - 2;
	             	Avy[valor1] = - Avy[valor1];
        	}
     	}
	/** si todos los cuerpos no estan emparejados hay que calcular los cuerpos sin parejas*/
	if(resto!=0){
		/** cogo el cuerpo en la posicion del primer cuerpo sin pareja mas mi nuemero de proceso si este pertenece al rango de los cuerpos sin pareja*/
		i=primero+yo;
		if(i<primero+resto){
          		ax[i] = fx[i] / m[i];
          		ay[i] = fy[i] / m[i];
	        	Avx[i] = vx[i] + ax[i] * interval;
	        	Avy[i] = vy[i] + ay[i] * interval;
	        	AUx[i] = x[i] + Avx[i] * interval;
        		AUy[i] = y[i] + Avy[i] * interval;
	        	if (AUx[i] <= 0) {
        			AUx[i] = 2;
	             		Avx[i] = - Avx[i];
        		}
	        	if (AUx[i] >= width ) {
        	    		 AUx[i] = width - 2;
	             		Avx[i] = - Avx[i];
        	  	}
	          	if (AUy[i] <= 0) {
        	  		AUy[i] = 2;
	             		Avy[i] = - Avy[i];
        	  	}
	          	if (AUy[i] >= height) {
        	     		AUy[i] = height - 2;
	             	Avy[i] = - Avy[i];
        	  	}
		}
		/** Si el resto es mayor que el numero de procesos aun habra cuerpos sin asignar*/
		if(resto>numeroProcesos){
			/** si la posicion dada por el primer porceso sin pareja mas el num de procesos mas mi id de proceso pertenece al rango coge el cuerpo
  				en esa posicion y lo calculo */
			if((primero+yo+numeroProcesos)<(resto+primero)){
				i=primero+yo+numeroProcesos;
				ax[i] = fx[i] / m[i];
        	  		ay[i] = fy[i] / m[i];
			        Avx[i] = vx[i] + ax[i] * interval;
		        	Avy[i] = vy[i] + ay[i] * interval;
		        	AUx[i] = x[i] + Avx[i] * interval;
		        	AUy[i] = y[i] + Avy[i] * interval;

			        if (AUx[i] <= 0) {
        				AUx[i] = 2;
	             			Avx[i] = - Avx[i];
	        		}
		        	if (AUx[i] >= width ) {
        		     		AUx[i] = width - 2;
	             			Avx[i] = - Avx[i];
	        		}
		        	if (AUy[i] <= 0) {
        		     		AUy[i] = 2;
	             			Avy[i] = - Avy[i];
	        		}
	          		if (AUy[i] >= height) {
        	     			AUy[i] = height - 2;
	             			Avy[i] = - Avy[i];
        	  		}
			}
		}
	}
	/** Marco el fin del tiempo de calculo, incremetno la variable total y marco el fin del evento en el MPE */
	tCalcF=MPI_Wtime();
	MPE_Log_event(stop_comp, 0, "end comp");
	tCalcT+=(tCalcF-tCalcSt);

	/** marca del inicio del tiempo de comunicaciones*/
	tComSt=MPI_Wtime();
	/** operacion de comunicacion de las variables locales al proceso raiz para pasarlas a la variable global*/
	MPI_Reduce(AUx,x,nbodies,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
	MPI_Reduce(AUy,y,nbodies,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
	MPI_Reduce(Avx,vx,nbodies,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
	MPI_Reduce(Avy,vy,nbodies,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);

	/** se comnunican las variables globales a todos los procesos*/
	MPI_Bcast(x,nbodies, MPI_DOUBLE,0,MPI_COMM_WORLD);
	MPI_Bcast(y,nbodies, MPI_DOUBLE,0,MPI_COMM_WORLD);
	MPI_Bcast(vx,nbodies, MPI_DOUBLE,0,MPI_COMM_WORLD);
	MPI_Bcast(vy,nbodies, MPI_DOUBLE,0,MPI_COMM_WORLD);
	tComF=MPI_Wtime();
	tComT+=(tComF-tComSt);

	/** se libera de la memoria las variables locales */
	free(ax);
	free(ay);
	free(AUx);
	free(AUy);
	free(Avx);
	free(Avy);


     return 1;
}
int main (int argc, char ** argv) {
   int i,iter;
     int x11 = 0;
     struct timeval start, end;
     double usec = 0;
     
     iter=0;
     if (argc < 4) {
        printf ("Usage:\n");
	printf ("   %s <input_file> <output_file> <X11_0|1>\n", argv[0]);
	printf ("Where:\n");
	printf ("   <X11_0|1> can be 0 or 1. 0 means no X11, 1 means use X11.\n");
        return 0;
     }

	/* Inicializacion de MPI*/	
     	MPI_Init(&argc, &argv);

	/* Inicializacion de MPE y de las varoables para mostrar tambie el tiempo de calculo*/
	MPE_Init_log();
	start_comp = MPE_Log_get_event_number();
	stop_comp = MPE_Log_get_event_number();
	MPE_Describe_state(start_comp, stop_comp, "COMP", "red");


	/* tomamos el valor del numero de procesos que usamos*/
	MPI_Comm_size(MPI_COMM_WORLD, &numeroProcesos);
	/* el proceso toma su valor de Id*/
     	MPI_Comm_rank(MPI_COMM_WORLD, &yo);

     if (!strcmp(argv[3],"1")) x11 = 1;
     init_bodies(argv[1]);
     if (x11) {
         allocate_grid(width, height);
         setupWindow(width, height);
     }

     for (iter = 0;iter < nsteps;iter++) {
         gettimeofday(&start, NULL); 
        if(strcmp(argv[4],"Cuerpos")==0){
	 	calc_nbodiesP();
	}else if(strcmp(argv[4],"Cargas")==0){
		calc_nbodiesC();
	}else if(strcmp(argv[4],"Secuencial")==0){
		calc_nbodies();

	}else{
		printf ("ERROR \n");

	}
         gettimeofday(&end, NULL);
         usec += end.tv_usec + 1000000*(end.tv_sec - start.tv_sec) - start.tv_usec;
         if (x11) {
             init_grid(width, height);
             for (i=0; i < nbodies; i++) {
                  grid[(int)x[i]][(int)y[i]]= m[i];
             }
             if(iter % REFRESH ==0 ) myDraw(width, height);
         }
     }


	if(yo==0){
     		if (x11) 
        		cleanup_grid(width, height);
	}
	if (yo == 0)
     		finalize_nbodies(argv[2]);

	/* Finalizamos el MPI */
     	MPI_Finalize();

	
	/* Imprimimos los tiempos de calculo, de comunicacion y totales de cada proceso*/
	printf("\n\n el tiempo de comunicacion para el proceso %i, es = %f \n ", yo, tComT);
	printf("el tiempo de calculo para el proceso %i, es = %f \n ", yo, tCalcT);
	printf("el tiempo total para el proceso %i, es = %f \n ", yo, (tCalcT+tComT));


     return 1;
}
